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1.0  INTRODUCTION 

The  purpose  of  this  investigation  was  to  improve  the  prediction  capability  of  Wright  Laboratory 
Three-Dimensional  Euler/N avier-Stokes  Aeromechanic  Method,  designated  TEAM,  through  the  imple¬ 
mentation  of  the  McDonnell  Douglas  nonequilibrium  Johnson-King  (J-K)  turbulence  model.  The  J-K 
turbulence  model  has  been  demonstrated  to  predict  more  accurately  surface  pressure  distributions  on 
supercritical  transonic  airfoils  and  wings  with  small  regions  of  axially  shock-separated  flow.  The  follow¬ 
ing  subsections  will  provide  a  brief  history  of  the  development,  calibration  and  limitations  of  the  J-K 
model. 

1.1  MDC  Johnson-King  Tbrbulence  Model 

The  McDonnell  Douglas  Research  Laboratory  has  developed  a  three-dimension  nonequilibrium  John¬ 
son-King  turbulence  model  in  support  of  the  ongoing  Delta  rocket  program  at  the  McDonnell  Douglas 
Space  Systems  (MDSSC)  Company,  Huntington  Beach,  California.  The  J-K  model  was  developed  to  pre¬ 
dict  more  accurately  the  surface  pressure  distributions  on  the  large  payload  variants  of  the  Delta  rocket. 
The  increased  payload  Delta  fbrebody  configuration  included  a  boat-tail  region  which  was  susceptible  to 
shock-wave-boundary-laycr  separation. 

Numerical  simulations  of  the  increased  payload  Delta  flowfield  with  the  standard  Baldwin-Lomax 
(B-L)  turbulence  model  underpredicted  the  extent  of  the  axial  separation  bubble,  lb  overcome  this  defi¬ 
ciency,  the  J-K  formulation  of  Abid  and  Johnson  [1]  was  implemented.  Simulations  were  performed  on 
tiie  Delta  N1B1  forebody  configuration,  shown  in  Figure  1.1.1,  with  the  Baldwin-Lomax  and  Johnson-King 
and  k-e  turbulence  models.  Figure  1.1.2  depicts  the  computed  surface  pressure  on  the  NIB  1  forebody. 
Note  the  improved  agreement  with  the  experimental  data  that  the  J-K  and  k-e  models  provide.  The  J-K 
formulation  developed  for  the  Delta  program  has  been  modified  for  implementation  into  the  TEAM  code. 
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Flgur#  1 .1 .2.  Comparison  of  Computed  and  Maasursd  Surfaca  Prassura  on  tha  Ml  B1  Forabody 
for  tha  MDC  Johnson-King,  Baldwin-Lomax  and  K-a  Turtouianca  Modal*;  M«=  0.8,  a  =  0.0, 

Rcs4,ooo,ooo.am. 


1.2  MDC  Johnson-King  Team  Code  Calibration  Plan 

To  demonstrate  the  J-K  implementation  into  the  TEAM  code  several  test  cases  were  computed.  The 
test  case  configurations  were:  (1)  a  simple  flat  plate,  (2)  the  RAE  2822  airfoil,  and  (3)  the  Onera  M6 
wing.  The  test  cases  were  chosen  to  provide  an  increasing  level  of  geometric  complexity.  For  all  test 
cases  the  analyses  were  conducted  with  both  the  B-L  and  J-K  turbulence  models.  Comparisons  were 
made  with  existing  experimental  data  where  possible. 

1.3  Turbulence  Model  Formulation  Limitations 

The  J-K  model  installed  into  the  TEAM  code  has  one  limitation  which  was  resulting  from  the  scope  of 
the  contract  as  spelled  out  in  the  statement-of-work.  The  J-K  model  in  the  TEAM  code  has  been  coded 
for  configurations  where  the  solid  walls  correspond  to  a  k=l  surface.  The  J-K  model  is  functional  for 
zonal  geometries;  however,  as  specified  above,  the  solid  wails  must  lie  on  a  K=1  surface. 
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2.0  GOVERNING  EQUATIONS 

In  the  following  sections  the  development  of  the  MDC  J-K  model  will  be  discussed. 

2.1  Johnson-King  Turbulence  Model  Formulation 

In  this  investigation,  the  J-K  formulation  of  Abid  and  Johnson  [I)  was  employed.  A  summary  of  their 
model  follows. 

In  the  Johnson-King  nonequilibrium  model  the  functional  form  of  the  eddy  viscosity  is: 


h  =  ho  (  1  ~  e-»  ) 


The  inner  viscosity,  pu.  is  defined  by  the  following: 


hi  =  Q  D2  x  N  x^ 


where  xm  is  the  maximum  Reynolds  shear  stress,  N  is  the  body  normal  distance  and  x  is  the  von  Karman 
constant  (x  =  0.4).  The  damping  coefficient  D  is  defined  by  the  following  expression: 


D  -  (  1  -  e&) 


where  A+  =  17  and  y+  is  the  standard  law-of-the-wall  coordinate.  The  Reynolds  shear  stress  is  assumed 
to  be  of  the  form: 


;  ;  h^ 


where  u>  is  the  magnitude  of  the  vorticity,  and  the  variable  xm  is  determined  from  the  solution  of  a  partial 
differential  equation  which  will  be  described  below. 


The  outer  eddy  viscosity  pt0  is  found  from  the  equation  below: 


ho  — ■  Q  <7  K  Cc  Fw  Yk 


where: 
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K-  0.0168 


where  6  is  the  boundary  layer  thickness  and  assumed  to  be  1.9  N^.  xmieq  is  the  value  of  assuming 
that  eddy  viscosity  profile  is  in  equilibrium  (i.e„  o  =  1.0 ).  The  diffusion  term  Dm,  modeled  by  Johnson 
[3],  is  defined  as: 
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Dm  - 


CD  A  I  1  - 


at(  0.7  6  -  Nm  ) 

Throughout  this  investigation  the  values  of  ai  and  Cd  were  fixed  at  0.25  and  0.5,  respectively. 
In  order  to  simplify  the  partial  differential  equation  a  change  of  variables  was  made.  Let: 


g  =  ^ 

=  4, 


,eq 


The  transformed  PDE  becomes: 


—  +  u  —  4-v  —  +  w  —  +r 
at  +  Um  ax  +  Vm  ay  +  Wm  az  +  1 


=  o 


where: 


(10) 


(11) 


(12) 


r  = 


al 

(A.  i 

\  Cd 

Lm 

1 

-4  ll 
-  a* 

2  Lm 

Igeq  1 

'  ai 

6  ( 

0.7 

-*)  J 

(13) 


The  above  equation  is  then  transformed  to  the  computational  plane  to  simplify  the  integration.  A  further 
simplification  is  made  by  assuming  steady  flow.  The  simplified  transformed  equation  becomes: 

[  um  sx  +  vm  £y  +  wm  ez]  ||  + 

!  «m  +  Vm  £y  +  Wm  +  i  =0 

*•  (14) 

Note  that  e  and  £  directions  correspond  to  the  transformed  streamwise  and  spanwise  directions,  respec¬ 
tively.  The  terms  containing  the  transformed  normal  derivative  ( r) )  were  assumed  to  be  negligible  along 
the  surface  where  xm  occurs.  For  complete  details  see  reference  [4].  The  above  equation  can  now  be 
solved  for  g  and  subsequently  xm.  Once  xm  is  known,  a  new  value  of  sigma  is  computed  from: 


^n+l  _  Qrnax  'tm 


(Ptlwl) 


max 


(15) 
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2.2  Solution  Procedure 

The  partial  differential  equation  for  g  (14)  was  integrated  via  successive-line-over-relaxation.  Line 
relaxation  was  performed  in  the  transformed  streamwise  (£)  direction.  The  cross-flow  derivative  (0 
terms  were  treated  as  source  terms  and  moved  to  the  right-hand  side  of  the  equation.  By  solving  the 
equation  tor  g  in  this  manner  overall  CPU  time  expenditures  were  reduced. 

When  the  TEAM  code  is  initialized  the  values  of  o  are  set  to  1.0.  Once  the  solution  has  been 
advanced  in  time,  the  JKTM  subroutine  writes  out  the  current  values  of  o  and  g  to  I/O  unit  55  (i.e., 
fort.55).  This  data  file  is  required  for  subsequent  restarts.  If  the  fort.55  data  file  is  not  present,  the  JKTM 
subroutine  defaults  the  value  of  o  back  to  1 .0. 


k 


i 
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3.0  MDC  MODIFICATIONS  TO  TEAM 

To  complete  the  installation  of  the  J-K  module  into  the  TEAM  code,  four  subroutines  were  added  to 
the  original  TEAM  software.  Subroutines  JKEXEC,  JKSET,  JKTM  and  THOMAS  are  discussed  below. 
Also,  a  minor  modification  was  made  in  subroutine  FILTER  to  allow  the  J-K  model  to  be  called  from  the 
main  program  driver. 

To  invoke  the  J-K  turbulence  model  the  turbulence  model  selection  term,  in  the  auxiliary  data  set,  is 
set  from  “bltm”  to  “jktm.”  Once  the  selection  term  is  set  to  “jktra,”  the  TEAM  code  will  invoke  the  J-K 
turbulence  model. 

3.1  Subroutine  JKEXEC 

Subroutine  JKEXEC  prepares  the  data  set  to  be  fed  into  the  J-K  turbulence  model.  The  function  of 
JKEXEC  is  to  search  through  the  boundary  condition  table  to  locate  the  viscous  wall  and  pass  the  vari¬ 
ables  to  subroutine  JKSET. 

3.2  Subroutine  JKSET 

Subroutine  JKSET  converts  the  flow  and  grid  data  passed  from  JKEXEC  into  the  form  required  by  the 
subroutine  JKTM.  JKSET  converts  the  cell-centered  to  cell-face  flow  variables.  Once  the  data  conver¬ 
sion  is  complete,  subroutine  JKSET  then  calls  subroutine  JKTM. 

3.3  Subroutine  JKTM 

Subroutine  JKTM  computes  the  turbulent  eddy  viscosity  in  a  manner  as  described  in  Section  2.0. 

3.4  Subroutine  THOMAS 

Subroutine  Thomas  performs  a  tridiagonal  inversion  required  by  subroutine  JKTM. 
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4.0  TEAM  CODE  CALIBRATION  TEST  CASES 

To  meet  the  contractual  obligation  of  the  statement  of  work,  MDC  performed  a  series  of  code  calibra¬ 
tion  runs  to  demonstrate  the  installation  of  the  J-K  turbulence  model  into  the  TEAM  code.  The  computed 
test  cases  were: 

la  --  Flat  plate  with  B-L  turbulence  model, 
lb  -  Flat  plate  with  J-K  turbulence  model, 
la  -  RAE  2822  Airfoil  with  B-L  turbulence  model, 

2b  --  RAE  2822  Airfoil  with  J-K  turbulence  model, 

3a  -  Onera  M6  wing  with  B-L  turbulence  model,  and 
3  b  --  Onera  M6  wing  with  J-K  turbulence  model. 

4.1  Flat  Plate 

Ti  iii  cases  la  and  b  were  computed  to  demonstrate  the  implementation  of  the  J-K  turbulence  model 
into  the  TEAM  code.  The  ffeestream  conditions  employed  for  this  test  case  were: 

M*  =0.7  Re  =  1,000,000 

Thi;  computational  grid  consisted  of  101  streamwise  by  61  normal  by  3  spanwise  grid  points.  A  com¬ 
parison  of  the  theoretical  and  computed  skin-friction  coefficient  for  both  the  B-L  and  J-K  turbulence 
models  is  shown  in  Figure  4.1.1.  Note  the  generally  poor  agreement  with  the  empirical  result  of  van-Dri- 
est.  Instead  of  decreasing  monotonically  down  the  plate,  the  values  of  Cf  increase  near  the  exit  of  the 
computational  domain.  This  behavior  indicated  that  the  TEAM  code  did  not  have  the  capability  to  accu¬ 
rately  model  near-field  subsonic  outflow. 


(Required  by  Contract) 
(Required  by  Contract) 
(MDC  Optional) 

(MDC  Optional) 
(Required  by  Contract) 
(Required  by  Contract) 
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Figure  4.1.1.  Comparison  of  computed  and  theoretical  skin-friction  coefficient  for  a  flat  plate  with 
the  Baldwin-Lomax  and  Johnson-King  turbulence  Models;  M  »  s  0.7,  Re  s  1,000,000.0. 

4.2  RAE  2822  Airfoil 

Test  cases  2a  and  2b  demonstrated  the  capability  of  the  J-K  model  to  compute  two-dimensional 
transonic  flow  over  a  supercritical  airfoil  section.  The  computational  grid  consisted  of  241  streamwise  by 
65  normal  by  2  spanwise  grid  points.  The  freestream  conditions  for  this  test  case  are  shown  below: 

M<»  =0.75  a  =  2.8  Degrees 

Rec  =  6.2X106 

Grid  point  spacing  in  the  normal  direction  was  maintained  such  that  the  y+  value  for  the  first  point  off  the 
wall  did  not  exceed  4.0.  The  TEAM  code  default  options  were  used  for  time  step  and  dissipation. 

Comparisons  of  computed  and  measured  surface  pressure  coefficients  are  made,  with  good  agreement, 
in  Figure  4.2.1.  The  J-K  solution  more  accurately  predicts  the  shock-wave  location  on  the  suction-surface 
of  the  airfoil.  Both  the  B-L  and  J-K  solutions  were  run  for  approximately  2500  iterations  and  4  orders  of 
magnitude  reduction  of  the  L2  norm  was  obtained. 
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Figure  4.2.1.  Comparison  of  the  Computed  and  Measured  Surface  Pressure  Coefficient  on  the 
RAE  2282  Airfoil;  M=°  =  0.75,  a  =  2.8  Deg.,  Re  =  6,200,000. 


4.3  Onera  M6  Wing 

Test  cases  3a  and  b  demonstrated  the  capability  of  the  J*K  model  to  compute  three-dimensional 
transonic  flow  over  wing  configurations.  The  computational  grid  for  the  ONERA  M6  wing  consisted  of 
193  stream  wise  by  36  normal  by  37  spanwise  grid  points.  The  freestream  conditions  for  this  test  case  are 
shown  below: 


M  «>  =  0.84  a  =  3.06  Degrees 


Rec=  11.7X106 

Grid  point  spacing  in  the  normal  direction  was  maintained  such  that  the  y+  value  for  the  first  point  off  the 
wall  did  not  exceed  10.0. 

Comparisons  of  computed  and  measured  surface  pressure  coefficients  are  made,  at  the  spanwise  loca¬ 
tions  of  (z/b)  =  0.2, 0.45, 0.65, 0.90  and  0.95  in  Figures  4.3. 1,2, 3, 4  and  5,  respectively.  The  agreement  of 
the  B-L  and  J-K.  computed  solutions  is  good;  however,  the  comparison  with  the  experimental  data  is  poor. 
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Figure  4.3.2.  Comparison  of  the  Computed  and  Measured  Surface  Pressure  Coefficient  on  the 
Onera  M6  Wing  at  the  y/b  =  0.45  Span  Station;  M*>  =  0.84,  a  =  3.0  Deg.,  Re  =  11,700,000.0. 
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Figure  4.3.3.  Comparison  of  ths  Computed  and  Measured  Surface  Pressure  Coefficient  on  the 
Onera  M6  Wing  at  the  y/b  =  0.65  Span  Station;  M  <*>  =  0.84,  a  =  3.0  Deg.,  Re*  11,700,000.0. 


Figure  4.3.4.  Comparison  of  the  Computed  and  Measured  Surface  Pressure  Coefficient  on  the 
Onera  M6  Wing  at  the  y/b  -  0.90  Span  Station;  M*  =  o.84,  a  =  3.0  Deg.,  Re  =  11,700,000.0. 


4-5 


-1.25 


Figure  4.3.5.  Comparison  of  ths  Computed  and  Measured  Surface  Pressure  Coefficient  on  the 
Onera  M6  Wing  at  the  y/b  s  0.95  Span  Station;  M«s  0.84,  a  =  3.0  Deg.,  Re  « 11,700,000.0. 


Poor  agreement  with  the  experimental  data  was  attributed  to  the  fact  that  the  computational  grid  did 
not  have  the  sufficient  density,  in  the  streamwise  direction,  to  capture  the  shock-wave  adequately  on  the 
suction  surface.  Grid  density  studies  were  required  but  not  possible  due  to  the  scope  of  the  contract  state- 
tnent-of-work. 

All  input  data  files  required  for  the  ONERA  Wing  simulation  are  listed  the  Appendix. 


4.4  Discussion 

MDC  has  implemented  a  J-K  model  into  the  TEAM  code.  Computations  have  been  performed  for 
both  2-D  fiat  plate,  2-D  airfoil,  and  3-D  wing  configurations.  Computed  solutions  from  the  B-L  and  J-K 
turbulence  models  appear  to  be  in  reasonable  agreement;  however,  for  the  case  of  the  ONERA  M6  wing, 
the  agreement  between  the  computed  solutions  and  the  experimental  data  is  poor.  Poor  agreement  with 
the  ONERA  data  is  attributed  to  inadequate  grid  resolution  on  the  suction  side  of  the  airfoil  in  the  recom¬ 
pression  region.  Further  calibration  of  the  TEAM  code  is  indicated,  but  not  covered  in  the  scope  of  this 
investigation. 
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6.0  APPENDICES 


Main  Program  Control  Input  Data  Set 
Auxiliary  Input  Data  Set 

‘zone  1’  ‘ _ tm’  0 

‘zone  T  ‘segment’  1 
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'zone  I’ 
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* 
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‘zone  1’  ‘segment’  1 
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24 
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l 

1 

Boundary  Condition  Data  Set 

‘zone  1’ 
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‘tlns-k’  9 

1  1 

1 
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33 

•far’ 

193193 

1 

37 

i 

33 

*far’ 

29  165 

1 

25 

l 

1 

‘solid’ 

1  29 

1 

25 

l 

1 

‘wake’ 

193165 

1 

25 

l 

1 

‘zone  1' 

165193 

1 

25 

i 

1 

‘wake’ 

29  1 

1 

25 

l 

1 

‘zone  1’ 

1  193 

25 

37 

1 

1 

‘wake’ 

193  1 

25 

37 

l 

1 

‘zone  1’ 

1  193 

1 

1 

1 

33 

'symmetry  _xz’ 

1  193 

37 

37 

l 

33 

‘far’ 

1  193 

1 

37 

33 

33 
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